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Variational principle of counting statistics in master equations 
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We study counting statistics of number of transitions in a stochastic process. For mesoscopic 
systems, a path integral formulation for the counting statistics has already been derived. We here 
show that it is also possible to derive the similar path integral formulation without the assumption 
of mesoscopic systems. It has been clarified that the saddle point method for the path integral is not 
an approximation, but a valid procedure in the present derivation. Hence, a variational principle in 
the counting statistics is naturally derived. In order to obtain the variational principle, we employ 
many independent replicas of the same system. In addition, the Euler-Maclaurin formula is used in 
order to connect the discrete and continuous properties of the system. 
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In a stochastic system, it would be possible to count 
a number of a certain target transition, and the num- 
ber of the target transition is a random variable. Since 
the random variable is directly related to a net current 
or flow in the system, the counting of number of transi- 
tions plays an important role in nonequilibrium physics, 
Especially, nonequilibrium properties for time-dependent 
systems have been studied largely; e.g., Brownian motors 
(ratchet syste ms) [H, [jl , H and a pump current problem 

6 a i ail Em 

There are several methods to calculate av erag e current 
in periodically time-dependent systems [l^, fislfl^l . The 
calculation of the average current enables us to have the 
no-pumping theorem and the pumping restriction theo- 
rem, which will become useful for the study of biochem- 
ical reactions or molecular machines [H, [3]. In addi- 
tion, recent progress gives recipes to evaluate not only 
the average current, but also all statistics including the 
fluctuation, i.e., 'counting statistics' 0, B B [HI ■ We 
note that basic concepts of full counting statistics (PCS) 
in condensed matter physics [TtI . [TsI . [19| are similar to 
the classical stochastic cases. 

For the counting statistics in master equations, there 
are mainly two different approaches. One is the tran- 
sition matrix approach [20l [2ll . [2^ , and the extensions 
of this approach have succeeded in calculating the pump 
current analytically 0, B, [13] ■ The other one is based 
on a path integral formulation Q . The transition matrix 
approach requires diagonalization of an effective evolu- 
tion Hamiltonian, which may be a complicated task for a 
system with a large transition matrix such as the ratchet 
problem. In contrast, the path integral formulation is 
available for the ratchet problem 

In [^, a path integral formulation for the counting 
statistics in mesoscopic systems has been investigated. 
The word 'mesoscopic' means that the system contains 
many particles or elements, so that the saddle point 



method in the analysis is approximately justifled. Hence, 
one can finally obtain saddle point equations similar to 
the Hamilton's canonical equations in classical mechan- 
ics. In the PCS, the similar discussions have been used 
[13, [il]- For the mesoscopic systems, a replacement of a 
discrete variable with a continuous variable (or replace- 
ment of a summation with an integral) may be justified 
because the discrete variable takes a large value. How- 
ever, if one loses the assumption of mesoscopic systems, 
the validity of the usage of the saddle point method and 
the replacement of the discrete variable with a continuous 
one are ambiguous. 

In the present paper, we will clarify that the saddle 
point equations similar to the Hamilton's canonical equa- 
tions are derived without the assumption of mesoscopic 
systems. Purthermore, it is clarified that the saddle point 
method is not an approximation, but a valid procedure 
in our derivation. Hence, we can finally obtain a kind 
of variational principle in the counting statistics. In or- 
der to obtain it, many replicas of microscopic systems 
are adopted, and the Euler-Maclaurin formula is applied 
in order to replace a discrete variable with a continu- 
ous one. This variational principle would become a basis 
for the counting problem, because it will allow us to use 
many analytical techniques developed in the Hamilton 
systems in order to investigate the counting statistics in 
stochastic systems. To demonstrate the derivation, we 
consider a simple model with 2x2 transition matrix in 
the present paper. Complicated models such as a ratchet 
problem could also be treated in a similar manner, by 
using a similar discussion in Q. 

Let us consider a simple model introduced in [7j for 
a pump current problem. In p], the similar system in 
the mesoscopic regime (the number of particles is very 
large) has been studied, but we here investigate the orig- 
inal system in which does not have the mesoscopic 
property. We here treat the following model: 



[L] ^ [bin] [R], 



(1) 
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where [L] and [R] are particle baths, and [bin] is assumed 
to contain either zero or one particles in it. The full 
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kinetic scheme is 

(i) [L] [bin] with rate ki{t); 

(ii) [L] ^ [bin] with rate 

(iii) [bin] — * [R] with rate k2{t); 

(iv) [bin] <— [R]with rate fc_2(i)- 

The transition rates {ki{t)} can be time-dependent, and 
in what follows, we denote them as {fci} for simplicity. 
The master equation is written as follows: 



reaction i in the iVmax systems is therefore 



d 
di 



p 




-fcl - fc_2 fc-1 + fc2 




P 






fcl + fc_2 — fc-1 - fc2 







(2) 



where p is the probability that [bin] is empty. Our goal 
is to evaluate the [bin] — > [R] flux for a finite time T. 

In order to calculate the counting statistics, we use the 
path integral technique. At first, we discretize the time 
a.st — mSt, m £ N. Hence, the final time is characterized 
by an integer M via T = MSt. The probability p is time- 
dependent, so that we denote the probability p at time 
t — mSt as pm in order to show the time dependence of 
p explicitly. 

We here define the number of reaction (or hopping) 
i G {1, —1, 2, —2} at time step m as AQm . A stochastic 
process described by a master equation can be simulated 
exactly by the Gillespie algorithm 2^ , and the num- 
ber of reaction i during time St obeys the Poissonian with 

average kiPm-i6t. For example, the probability of AQm' 
may be given by 



P(AQW|p„_i) = exp(-fcipm_i(5i) 



(kiPm-iSt) 



{AQ 



(3) 



An occurrence of a reaction changes the state, but the 
state variable p is continuous and AQm is discrete. 
In order to connect these quantities, we discretize the 
continuous variable Pm G [0, 1] as p,„ = rimSn, where 
n„i G {0,l,2,...,7Vmax} and 5n = l/A'^max- Note that 
we can recover the continuous property of the probabil- 
ity Pm when we take 5n ^ (i-e., -/Vmax oo). This is 
the most important point in the present paper; we should 
finally take iVmax — cxd in order to recover the continuous 
property of pm- 

There is still one problem as follows. In the Gillespie 
algorithm, one reaction changes the state of the system, 
and then the next reaction should be evaluated by using 
the updated state. Hence, only one reaction step, i.e., 
AQm = or 1, should be allowed. However, the Poisson 
statistics (eq. (jS])) does not include such effects. In other 
words, we should restrict AQm to {0, 1} for the natural 
interpretation of the Gillespie algorithm, but the domain 
of ([3]) is {0, 1, . . . } = N; it seems inadequate to use ([3]) 
directly. 

In order to avoid this problem, we introduce many in- 
dependent replicas of the system with the same state 
Pm-i (or rim-i) at each time step; totally A^max systems 
with the same initial state rim-i- The total number of 



AQ« ^ £ AQ«, 

and the probability of AQ^m is given by 

(kiUm-iSty 



(4) 



P(AQ^Vm-i) = exp(-fcinm-i(5t) 



(AqW)! 



(5) 



Here, we used the independent property of each replica. 
{AQm} for i = —1,2,-2 are obtained by the similar 
manner. Note that AQm' is the total number of reac- 
tion 1 in A'^niax replicas, and therefore AQm' G N (not 
{0, 1}). If we choose 6t small enough, it would be possi- 
ble to assume that there is at most one reaction in each 
replica (after taking A^max ~* c»). The introduction of 
the independent replicas, therefore, enables us to avoid 
the above problem (i.e., AQm should be or 1). In addi- 
tion, we will see later that the introduction of the replicas 
is necessary to use the saddle point method. 

In order to count the flux [bin] — > [R] , we introduce the 
quantities Q^ = AQ^' - AQ^'' and = AQifJ - 

AqI^^' , and calculate the characteristic function of Q^^ 
during the time interval T. The characteristic function 
conditioned by the initial state no is calculated from 

E[e*xc(Qf+Q?+---+Q5f-i+QM)|^(j] 
= E[e*xcQ?E[e*(xcQ?+-+xcQ5.)|j,^]|„g] 



= E[e*^<=^?E[e*'^°'2^E[- • • E[e*xcQ^j | . . 



(6) 



where xc may be assumed to be real. In general, a char- 
acteristic function is rewritten in a form of exp[S'(xc)] 
(see (|14l) and ([15])). When we define a quantity Q^ = 
J2m=i Qnn derivatives of S{xg) give cumulants for Q^; 
e.g., (Q^) = -idSixc)/dxc\xc=o, (5'Q^) ^ ((Q^)') - 
(qR)2 = {-i)^d^Sixc)/dx\c=o, etc. The quantity 
corresponds to the fiux [bin] [R] in the A^max systems, 
and hence an explicit calculation for the characteristic 
fimction ^ is needed to obtain the counting statistics 
for the fiux [bin] ^ [R]. 

For simplicity, we here show a calculation for 
EjgixcQm ^}jg calculation of the characteristic 
function ^ can be performed in a similar manner. By 
using the Fourier transformation of the Kronecker delta. 



S{A,B) 



dx 
2^' 



(7) 
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we obtain 
E 



-| N^max OO OO OO OO 

E E E E 



in ([8|). Hence, the second term in ((TO)) vanishes when 
-^max ~* 00, and we can verify the replacement of the 
discrete variable with the continuous variable. Note 
that the domain of the integration becomes [0, 1] when 
■^max — *■ oo; this is consistent with the fact that p is the 
probability. 

Finally, we obtain the characteristic function ^ in the 
path integral form as 



n„,— 



exp [kinyn-ie+xrn.^A 



X exp [k2 (-/V, 
X exp [fc_i(7V, 

X exp [fc_2 



max ^m — 1 — Xtti+XC 
max '^m—l)^ — Xm^'^ 



(8) 



where e±i 



1. Note that the probabilities 



P(A(5m|nm-i) for « = ±1,±2 are correlated via the 
Kronecker delta, and then a normalization constant Zm 
is introduced. 

Equation ([8]) has still one summation for rim, and its 
discrete property is inconvenient for further analytical 
treatments. However, different from the mesoscopic case 
in it may be inadequate in the present case to sim- 
ply replace the summation for rim with an integral. We 
therefore use the Euler-Maclaurin formula for a function 



1 "2 
Tl— Til 

r(n2 + l/2)/L 



g{x)dx 



(ni-l/2)/L 



1 



24L2 



, /n2 + l/2\ , All -1/2 



L 



L 



(9) 



where g'{x) is the derivative of with x. In our calcu- 
lation, we have the following form with a function f{n): 



l + l/(2Af„,x) 



Tl = 

1 

~ 24 



/'(I 
1 



l/(2W„ax) 
1 



2A^n 



2A^n 



=,Wn,ax/(- 



(10) 



where f{x) is determined by ([5]) (or eventually, by ©). 
Although f{x) is a complex function, we can easily see 
that the real part of f{x) is negative because Re e±^ < 



n / , """^ '^p- 



2tt 



,AfmaxS(xC,{p™}:{X™}) 



(11) 



where 

S{XcAPrn},{X7n}) 
M 



m— 1 



Pm Pm — 1 N 
H{xC,Pm~l,Xrn) 



(12) 



and 



^^(XC,Pm-l,Xm) 

= - ^lPm-ie+x„ - ^2(1 - Prn-l)e-x,r.+XC 



(13) 



Note that A^^ax — > cx) is a necessary condition in order 
to recover the continuous property of {pm}- Hence, we 
use the saddle point method and obtain 



E[e'XC E„ Ql Ipq] ^ exp [A,„ax5(XC, {PmV, {XmV 



(14) 



In addition, the random variables {Q^} independent 
each other, so that 



Pol 



exp[S{xcAP^nr\{Xmr')] , (15) 



where the superscript 'cl' means that these quantities are 
determined by saddle point equations (see (fT7|) and (fT8|) ). 
After taking the continuous time limit, we finally obtain 
the following result 



E[e'xc/d"dtQ^(t)|p^ 
~ exp 



H{xc,p'Ht),x'\t)) 



(16) 



where x'^'(i) and p^^{t) are evaluated by the saddle point 
equation 



.dp'Ht) _ , dHixc,p'\t),x''\t)) 



dt 



dx-\t) 



.dx'\t) _ dH{xc,p'\t),x'\t)) 



dt 



dp^^t) 



(17) 
(18) 
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Different from the mesoscopic case in [8| , the saddle point 
method is not approximation, as explained above; (|16p . 
(fT7| and (fT8|) are valid because we finally take A^max 
cxD. As a result, we obtain a kind of variational principle 
in the counting statistics. Namely, when we define the 
action S as 

the path of x and p are given by the variational prin- 
ciple, and then we have equations similar to Hamilton's 
canonical equations ([17]) and (fT8|) . 

Using the variational principle, we can calculate a 
pump current when the transition rates ki and k_2 vary 
adiabatically with a period Tq (fc_i and k2 are assumed 
to be time-independent). Following the similar analyti- 
cal treatment in [8], we define Sc as the surface enclosed 
by the contour c in the space of parameters ki and fc_2, 
and set the time derivatives in ([T7]) and (fT5|) as zero due 
to the adiabatic condition. We finally obtain 

E[e*xc/<rdtQ-(t)|p^] _ gS,e„„,+5.,^ (20) 

where 

Sgcom^-^<f dfcl<ifc-2-Ffci,fe_2, (21) 



Here, we define = ki + fc_2 + ^-i + ^2 and K = 
(fc^4fcifc2e;^c + 4fc_i/c-2e-xc)^^^- This is consistent with 
the results in Q. 



In conclusion, we derived a variational principle for 
a counting statistics in a microscopic master equation. 
The word 'microscopic' means that we do not need an 
assumption of the mesoscopic systems. Due to the lack 
of the mesoscopic properties, we cannot simply connect 
discrete variables and continuous variables. Hence, the 
derivation of the path integral formulation was carefully 
performed. In the derivation, the continuous state vari- 
able is replaced with a discrete variable, and the con- 
dition to recover the continuous property validates the 
usage of the saddle point method. One of advantages of 
this variational scheme is that we may use many analyt- 
ical method in Hamilton systems in order to study the 
counting statistics in stochastic processes. While further 
studies will be needed in order to make the variational 
principle given in the present paper useful, the present 
work will give a basis for future works of counting statis- 
tics. 



Fk,,k^, = -e-xc (e'^°fc2 + (22) 

and 

/"To 

Sci^^ dt{K+ - K). (23) 
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